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ABSTRACT 

We present a Bayesian Voronoi image reconstruction technique (VIR) for in- 
terferometric data. Bayesian analysis applied to the inverse problem allows us to 
derive the a-posteriori probability of a novel parameterization of interferometric 
images. We use a variable Voronoi diagram as our model in place of the usual 
fixed pixel grid. A quantization of the intensity field allows us to calculate the 
likelihood function and a-priori probabilities. The Voronoi image is optimized 
including the number of polygons as free parameters. We apply our algorithm 
to deconvolve simulated interferometric data. Residuals, restored images and 
values are used to compare our reconstructions with fixed grid models. VIR has 
the advantage of modeling the image with few parameters, obtaining a better 
image from a Bayesian point of view. 

Subject headings: methods: data analysis — methods: numerical — methods: 
statistical — techniques: image processing — techniques: interferometric 



1. Introduction 



Astronomical interferometric data result from the addition of instrumental noise to the 
convolution of the sky image and the instrumental response. Because of incomplete sampling 
in the {u, v) plane, obtaining sky images from interferometric data is an instance of the inverse 
problem, and involves reconstruction algorithms. 
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The CLEAN method consist s of modeling t he side-lobe disturbances and subtract them 
iteratively from the dirty map (iHogbonj 119741 ) . The CLEAN method works well for low 
noise and simple sources. But if the source has many complex features , or if the data is too 
noisy, CLEAN will do only a few iterations returning a noisy image (IHogbonj Il974j ). An- 
other shortcoming is that CLEAN involves some ad-hoc parameters (the loop gain, stopping 
criteria, clean beam) that bias the final reconstruction, in the sense that CLEAN can give 
many different reconstructions for the same dataset. 

The maximum entropy method (MEM) finds the image that simultaneously best fits 
the data, within the noise level, and maximizes the entropy S. This is done by minimizing 
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where, for the case of interferometric data, can be calculated as 
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where the sum runs over all the AVis visibilities, the symbol ||2;|| stands for the modulus 
of the complex number z and is the root mean square (rms) noise of the corresponding 
visibility. A is a control param eter and the entropy S varies for different implementations (e.g. 



M. 



Narayan fc Nityanandal 119861 ). The entropy is used as a regularizing term in a degenerate 
inverse problem, when there are more free parameters than data. Different formulations for 
S appear in the literature. Some examples are J2i Si lii(-^i)) y^,jln(pj,), p,; Infpj) 

where Jj is the specific intensity value at pixel i and Pi 
and references therein). 
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Cornwell fc Evand (1l985l ) used MEM in the AIPS VM task. Their method makes some 



approximations that diagonalize the Hessian matrix required to optimize their merit function. 
They used an entropy of the form 5* = — /j log (/j/mj), where the sum extends over all 
the pixels i, is the model image and {mj}"^^ is a prior image. However, the neglect of 

the side-lobe contribution to the Hessian may lead the optimization to local minima that still 
bear instrumental artifacts. ICasassus et al.l (120061 ) implemented a MEM algorithm based on 
the conjugate gradient method, without the use of the Cornwell and Evans approximation. 
They used an entropy of the form S = — ^ ■ Jj log (li/M), where {Ii}^i is the model image 
and M is a small intensity value, i.e they start with a blank image prior, and M is an 
intensity value much smaller than the noise. 

Bayesian analysis is a powerful tool for image reconstruction techniques. In this applica- 
tion, our goal is to find the most probable image by maximizing its a-posteriori probability. 
For Bayesian methods, the a-priori and likelihood distributions are needed. To derive the 
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a-priori probability the definition of an intensity quantum is needed. Tliis quantum repre- 
sents tfie minimum measurable intensity unit. The intensity in each pixel can be interpreted 
as a number of quanta = a^Ni, where is the intensity in pixel i, Uq is the quantum size 
and Ni the number of quanta in pixel i. 



Pina fc Puetteii (119931 ) used Bayesian analysis in the Pixon algorithm. They use a 
variable model and maximize P{I, M\D), that is, the probability of the image / and model 
M given the data D. In their approach the model used to parameterize the image is a 
set of Gaussians which are used to average a pseudo-image. The pseudo-image starts as a 
maximum residual likelihood reconstruction and a local Gaussian pixon is assigned to each 
of its pixels. The number of pixons, and hence the number of free parameters, is reduced in 
each iteration. 



Sutton fc Wandeltl (120061 ) have used Bayesian analysis for interferometric data, but using 
a fixed pixel grid to parameterize the model image. They use Gibbs sampling to determine 
the posterior density distribution. 

The most typical model used in astronomy to represent the sky brightness distribution 
consists of a pixel grid. A big disadvantage of this grid is that the number of pixels remains 
fixed as well as their size. Often, uniform pixel grids involve more free parameters than really 
needed to fit the data. 

The purpose of this paper is to explore Bayesian reconstruction with image models based 
on Voronoi tessellations in place of the usual pixelated image. We call this new deconvolution 
method "Voronoi image reconstruction" (VIR, hereafter). The advantage of using Voronoi 
models is that it is possible to use a smaller number of free parameters, as required by 
Bayesian theory. Our purpose is not optimal CPU efficiency; we search for the optimal 
image and model from a Bayesian point of view. 



We used the Cosmic Background Imager (CBI, iPadin et al.l |2002| ) to illustrate our 
method. The CBI is a planar interferometer array with 13 antennas, each 0.9 m in di- 
ameter, mounted on a 6 m tracking platform. An example of CBI baselines is shown in 
Figure [H The radius of the hole at the center of the [u, v) plane is the reciprocal of the 
minimum distance between two antennas, measured in wavelengths. The side-lobes of the 
CBI are caused mainly by this central hole in the {u, v) baselines. 

We briefly summarize the elements of Bayesian theory that determine the probability 
distributions concerning our problem (Section [2]). The new model based on Voronoi tes- 
sellations is described (Section [3]), as well as optimization issues involved in our problem 
(Section H]). We discuss implementation details such as the optimal quantum size and num- 
ber of Voronoi polygons (Section [5]), compare reconstructions made with MEM and VIR 
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Fig. 1. — Coverage in the {u, v) plane of the CBI in the configuration used for our simulations. 
(Section [6]) and finally summarize our results (Section [7]). 

2. Bayesian Theory 

An image model is required to parameterize the sky brightness distribution. The most 
typical model used in astronomy is a rectangular grid of uniform pixels. That configuration 
of pixels is the model M, and the distribution of brightness in the model is called an image 
/. We search for the image that represents as accurately as possible the visibility data D. 
The Bayesian image reconstruction approach, using a fixed model, tries to find the image 
that maximizes the probability P{I\D, M), i.e. find the most probable image given the data 
and the model. 



Using the Bayes theorem, we obtain 



PiI\D,M) 



P{D\I,M)P{I\M) 
P{D\M) ' 



(3) 



Since the data is fixed, P{D\M) is a constant in the problem when the model is not considered 
as a variable. Thus, the fixed image model optimization problem reduces to 



maxP(/|D,M) = max P{D\I,M)P{I\M). 



(4) 
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The first term, P{D\I, M) is called the likehhood, and measures how well our data represents 
our image. The second term, P[I\M) is called the image prior, and gives the a-priori 
probability of the image given the model, i.e. how probable is the image given only the 
model. 

In the case of having a variable model, what we would like to find is the image and 
model that maximize P{I, M\D), i.e. find the most probable image and model given the 
data. In this case we find 

P{I,M\D) = P{I\D,M)P{M\D) 

P{D\I,M)P{I\M)P{M\D) 

P{D\M) 
P{D\I,M)P{I\M)P{M) 



P{D) 



(5) 



Since the data is fixed, P{D) is constant in our problem. As we cannot privilege one model 
over another in the absence of image and data, P{M) is the same for all models, so it is not 
important for our analysis. This way, our optimization problem reduces to 

maxP(/, M\D) = maxP(D\L M)P(I\M). (6) 



2.1. Probability Distributions 

Our data is a set of A^vis observed visibihties {V^^^, V^^^, • • • , V^^l}. If we have a certain 
model M and image /, we obtain model visibilities {V^™°'*} by simulating the interferometric 
observations over our image: 

A{x, y)I{x, y) exp [2m{ukX + v^y)] / , (7) 

-oo yl-x^-y^ 

where {uk-, Vk} are the coordinates of baseline k in the (m, v) plane and A is the primary beam. 
We thus have a set of AVis model visibilities. Assuming that each visibility is independent 
from the others and Gaussian noise, the likelihood is 



rV[iod\ 
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To obtain the image prior, P(/|M), we calculate the statistical we ight of a given dis- 



tribution of counts (as in lPina fc Puetterlll993l : ISutton &: Wandeltll2006l ). Consider a model 



consisting of n cells. In the case of a traditional image, each pixel would be a cell. There is 
a number of N quanta falling into these cells. These are intensity quanta of some size Uq. 
In the case of a pixelated image, the intensity in each pixel i would be U = a^Ni, where Jj is 
the intensity in cell i. Each quantum could fall into any of the n cells, so the total number 
of possible configuration for the N quanta will be . The probability of the image given 
the model is the probability of a certain state {A^^i, N2, ■ ■ ■ , A^„} that represents that image, 
where iVj is the number of quanta in cell i. Consider a given image configuration defined 
by a particular distribution {Ni}. The image distribution is not changed in the N\ possible 
redistributions of counts between cells, provided each Ni is constant. The Yli swaps of 
counts within each cell keep the same image configuration. The model M consists of the 
Voronoi diagram and the total number of quanta (i.e. n, the position of the generators and 
N), thus the a-priori probability is 



P{I\M) = P{{N,}\n,N) 



n 



N 



(9) 



As explained above, CTq is an intensity quantum. It is also possible to describe the 
number of quanta per cell using a flux quantum erf, where i is the index of the cell to which 
we associate the quantum. This flux quantum can be expressed in terms of the intensity 
quantum as af = (XqAj, where Ai is the area of cell i. In this case, the number of quanta 
per cell is Ni = Fi/af, where = liAi is the flux of cell i. This leads to Ni = h/a^, 
which is the same expression for A^j obtained using the intensity quantum (jq. Using these 
cell- dependent flux quanta, the probability of a quantum falling into each cell will be ^ for 
every cell, leaving the a-priori probability the same as Eq. [9l 



2.2. MEM and Natural Entropy 

In Bayesian theory, for a fixed model, the image / can be found by optimizing the 
a-posteriori probability: 

maxP(/|D,M) = mm{- \nP{D\I,M)P{I\M)) 

= mm > ^-^ In ^ 

= mm^x^-5', (10) 
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where we have defined the natural entropy S* = In i n^YlNi ) ■ ISutton fc Wandeltl (120061 ) call 



the term In (A^!/ Hi multiplicity prior. In the limit of large Ni. 

S ~ A^ln 2^^^^^^^' 1^^^) 

i 

and it can be seen that the Bayesian method is very similar to MEM in the sense that we 
are adjusting the image to the data while maximizing an entropy of the form of Eq. [Ill VIR 
uses the natural entropy as a regularizing term. 



3. A New Image Model based on Voronoi Diagrams 

A Voronoi diagram is a division of the Euclidian plane into n regions Vj defined by n 
points Xi (called sites or generators) such that every coordinate x in the space belongs to Vj 
if and only if — afjH < ||x — V j ^ i. The result of the above definition is a set of 
polygons defined by the generators. Fi gure [2] shows an exa mple of a Voronoi diagram. For 
further details on Voronoi diagrams see lOkabe et al.l (119921 ). 




Fig. 2. — Example of Voronoi diagram. 

We propose a 2D Voronoi diagram in place of the usual pixelated, uniform grid, image 
as our model. We associate an intensity Jj to each of these polygons. The advantage of using 
a Voronoi diagram is that we can use just as many cells (i.e. free parameters) as the data 
requires. Our optimization parameters will be the position of each generator Xi = {xi,yi), 
and the intensity at each cell, Jj. 

With our new model M consisting of n generators {3 x n parameters, Xi,yi and Jj for 
each generator), we can vary the number of free parameters as required by the optimization 
problem. We can see in equation (fTOl) that the entropy S increases as the number of cell n 
decreases. 
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4. Optimization 

The optimization problem can be seen as a maximization of the a-posteriori probabihty 
max/^M M\D), or equivalently as a minimization of the more convenient merit function 
L = — 5*. The conjugate gradient method (CG) is often used for minimization problems 
where derivatives can be easily calculated. Though it is usually fast in convergence, CG has 
the problem of converging on local minima depending on the initial condition. The use of 
other optimization algorithms is postponed to future work. 

The CG method searches parameters space using the gradient of the function to be 
minimized. The derivatives of this function are 

dx 2 dx dx' 

2 ^Vis -I / QT/mod 



dx 



dx 



2 E ((^r' - vtr^) , (13) 



where x is any of the optimization parameters {xi, t/i or Jj). The derivatives of the visibilities 
with respect to the position Xj = {xi,yi) of the i generator are 



mod 



dxi 



gyraod 



dyi 



J2 [ih - h) Yl ^i^tiiMJi + 6,)e(*'^^+^«^i)] , (14) 

j€Ji /|pixcl leaij 

E AM(M,t, + ye(*'^^+^°^^)], (15) 



jdiJi / 1 pixel lea 



where Jj is the intensity in cell i, J, is a set of the indices of the polygons adjacent to Vj, Ojj 
is the edge which divides polygons Vj and Vj, I sums over the pixels which intersect aij, A is 
the CBI primary beam. For further details see Sec. \M 

The derivative of the visibilities with respect to the intensity of each cell Jj is 

gyraod ^ ^-^ ^nUkAx) siu (nVkAy) y ^^^2.riu,cc,+v,n) ^ (IQ) 

dli Tc'^UkVk ^ ' 

pixels leVi 

where fc^ = (u^, is the baseline corresponding to the pair of antennas fc, Ax and Ay are 
the pixel width and height, and the sum extends over all the pixels inside Vj. 

The entropy only depends of the intensities Jj, so ^ = ^ = 0, then (see Sec. IA.2P 



fc=l k=\ 
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5. VIR Design and Implementation 

We have designed, and implemented in C++, VIR with 6 modules which include algo- 
rithms for: 

• the generation of the Voronoi diagram 

• calculation of model visibilities 

• calculation of the merit function L to be optimized as well as its derivatives 

• fitting a Voronoi diagram to an image 

• the CG method 

• the optimization of the number of polygons 



VIR uses the CG method from iPress et al.l (119921 ) and searches for the position and 



intensities of the Voronoi polygons, Xi,yi,Ii, that minimize our merit function L. The CG 
method modifies the intensities and also moves the positions of the Voronoi generators. 
This causes the shape of the Voronoi polygons to change as well. A general problem with 
CG is that it usually converges on local minima. For VIR in particular, though Voronoi 
polygons intensities adjust quite fine, the positions of the generators are difficult to modify 
substantially. The VIR parameter space is smooth enough in intensity space to converge to 
a good solution. But the parameter space in cell generator positions is very structured, and 
CG is quickly stuck on local minima. 

Due to the fact that CG easily falls into local minima, we needed a good approximation 
for the initial Voronoi diagram. For this purpose we used a pixelated version of the Bayesian 
algorithm, where the model was a uniform grid. We decided to do a pure (maximum 
likelihood, ML) reconstruction and use the fifth CG iteration as our starting image. We 
chose this particular iteration because on inspection the modeled images were still smooth. 
Pure reaches convergence with noisy images, where the true image is unrecognizable. We 
then fitted a Voronoi diagram to the image (see Sec. [B]) and ran CG using the positions and 
intensities of the generators as our free parameters, which led to our final reconstruction. 
Truncation to a level of 10~^ quanta was used to enforce positivity. 



An important issue to consider is the size of the quantum a^. ISutton fc Wandeltl (120061 ) 
treat free parameter. But, as we now explain, was held constant in this imple- 

mentation of VIR. We treat the number of quanta per cell as a continuous variable in order 
to use the CG method. Entropy is maximized at (jq = oo, where, for a given configuration of 
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intensities {/?,}, = and 5 = 0. For every other value of A^, the entropy will be negative. 
This means that even for large cXq, the intensities Jj = (J^Ni can have reasonable values (using 
small Ni). Figure [3] shows as a function of for 51 Voronoi generators and 3 different 
intensity distributions using the model tessellation of Figure SJ:;. We considered: 1- the VIR 
intensities of Figure Ht, 2- a uniform intensity distribution image [Ni = ^ V i), 3- a spike 
where all are only in one cell (A^, = N,Nj = 0\/j^i). The curves of Figure [3] are 
obtained by keeping the intensities fixed and modifying cXq in order to obtain different A^. 
It can be seen on Figure [3] that the entropy is maximized at A^ = 0, independently of the 
intensities {/»} of the model, where the optimal value of Uq = cxd is achieved if the number 
of quanta per cell is tre ated as a continuou s varia ble. If the number of quanta per cell were 
discrete variables, as in ISutton fc Wandeltl (120061 ) . the choice of a big ciq would admit only 
zero values for every cell. Otherwise, if one or more quanta fell in a given cell, the intensity 
of that cell would diverge as for arbitrarily large Cq, causing a big value. Therefore, in 
our continuous optimization the intensity quantum must be determined a-priori. 




Fig. 3. — Entropy values for different A^, n = 51 and keeping {Jj} fixed. This is achieved by 
varying cTq. (a) VIR reconstruction intensities, (b) Uniform intensities distribution, = ^ 
V i. (c) Only one cell has all the quanta. = N,Nj = 0\fj^ i. 

In the Bayesian description of the entropy we count events that fall in each cell. It seems 
reasonable to take the noise level as the minimum value of intensity we can distinguish. So, 
CTq should approximate the estimated therma l noise in the natu rally weighted dirty map. 



Briggs et al.lll999l ) is 



The definition of the weighted dirty map (e.g. 

/oo /"OO 
/ iy(M,t;)\/(M,tOe-2-*("^+''J')rfurft;, 
-oo J ~oo 

W{u,v) = — y^^Wk5{u -Uk,v - Vk), 



(18) 
(19) 
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Fig. 4. — Comparison of MEM and VIR reconstruction techniques for a SNR of ~ 52. (a) 
The true image, (b) Dirty map. (c) VIR reconstruction with its polygons drawn, (d) VIR 
reconstruction, (e) Dirty map of the VIR reconstruction residuals, (f ) Restored image for the 
VIR model, (g) MEM reconstruction, (h) Dirty map of the MEM reconstruction residuals, 
(i) Restored image for the MEM model. 



where the sums extend over all visibilities, are the weights given to visibility k and 5 
is the two-dimensional Dirac delta function. Propagating the thermal noise, we get for the 
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standard deviation of the dirty map 



rms 



D 



2' 



(20) 



where cxfc are the visibihties standard deviations. To take into account model pixels correlated 
by the interferometer beam, we should multiply the previous expression by ^/Nheam, where 
-^beam IS the number of pixels inside a beam pattern. This leads to 



We calculated the noise with natural weighting, Wk = because this is the weight we give 
to each individual visibility data in the optimization of the merit function. 

Once we have the value of cxq we search for the optimal number of cells n. In Figure [5] 
we plot the optimal merit function for different n and cXq. These reconstructions were made 
over a simulation of CBI observations on a mock sky image (Figure UK). We averaged over 
100 reconstructions with different realizations of Gaussian noise. The average curves shown 
in Figure [5l start with n = 10 and end with n = 100 for even n. One single reconstruction 
for all n took about two hours using an AMD Athlon64 XP3000 processor with 1GB of DDR 
RAM at 333 MHz, so the 300 reconstructions took about 25 CPU days, but we distributed 
the work in six computers, so it took about 5 real days in total. It can be seen that for a 
signal to noise ratio (SNR) of ~ 52, on average, the optimal number of polygons n is between 
50 and 55. When Cq is diminished to jqCt^, on average, the optimal merit function is found 
at n close to 30. For cTq = lOciq, the optimal n is found between 80 and 90. It can be seen 
that as we increase the value of we reach lower values for our function, as discussed above. 
Furthermore, the optimal number of polygons increases. 




(21) 



For natural weights, = 




(22) 



6. Example Reconstruction 



6.1. 



Mock Dataset 



The mock sky image we used for simulations is a 256 x 256 image consisting of three 
Gaussians and a rectangle. Figure H^i shows this image on a 128 x 128 pixel field. Pixels 
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Fig. 5. — The merit function L for different Uq and number of polygons n. The hues are 
averages taken over 100 different reahzations of noise for each n. (a) Reconstructions made 
using (Tq = ^^cTrms- (b) Recoustructious made using cr^ = cr^rns- (c) Reconstructions made 
using (Tq = lOdrms- (d) L as a function of n for a practical application of VIR to the simulated 
visibilities used in the reconstructions of Figure HI In this practical application, the minimum 
L was found at n = 51, and is indicated by a vertical line. 



are 0.75' x 0.75', while the CBl's primary beam is of 45' FWHM (60 pixels), so most of the 
emission lies under the beam. We simulated a CBl observation of 3620 visibilities over this 
image and added Gaussian noise to the visibilities in order to reach a SNR of ~ 52. This SNR 
was calculated by taking the maximum intensity from the dirty map using natural weights, 
and using the noise (see Eq. |20]). Simulation of the CBI observations is performed 
with the MockCBI program (Pearson 2000, private communication), which calculates the 
visibilities V{u, v) on the input images I{x, y) with the same uv sampling as a reference 
visibility dataset (Eq. [7]). Thus MockCBI creates the visibility dataset that would have 
been obtained had the sky emission followed the true image. Figure H b shows the dirt y map 
calculated over these simulated visibilities using the DIFMAP package (jShepherdll 1 99 71 ) . The 
CBI's primary beam is drawn as a dashed circle. The secondary side-lobes due to the central 
discontinuity in u-v coverage can be distinguished in Figure at a level comparable to the 
true emission. 



6.2. MEM Reconstruction 



The VIR method was compared with the MEM algorithm described in ICasassus et al. 

( 120061 ). To fit the model image to the observed visibilities, MEM calculates the model 
visibilities required by its merit function Lmem- The model visibilities are those obtained 
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by a simulation of CBI observations had the sky followed the model image . The free- 
parameters of our MEM model are the pixels in the model 64 x 64 image. The model 
functional we minimize is Lmem = ~ ^S, with the entropy S = — ^-/jlog/j/M, where 
M is a default pixel value well below the noise level, and {/j}^^ is the model image. We 
started with the fifth iteration of a pure reconstruction (A = 0) as initial condition for the 
CG minimization. This is the same ML initial condition used in our VIR method. Figure 
shows the reconstructed image using A = ^ and M = IQ-V^ms inset on a larger 128 x 128 

0" rms 
_ 

6.3. VIR Reconstruction 

The MEM algorithm described above requires the prior assignment of the A and M 
parameters as well as the entropy formula. In contrast, our VIR algorithm is free from such 
arbitrary parameters (provided the optimal cxq is indeed equal to o"rms)- For our VIR method, 
we only need to find the number of polygons to be used. In order to find the optimal number 
of polygons we reconstructed with different numbers of generators in a range covering each 
natural number from n = Q to n = 100. We found a minimum at n = 51. Figure [5] 
summarizes this search. The whole search for a particular realization of noise took about 10 
hours on the AMD Athlon64 XP3000 processor with 1GB of DDR RAM at 333 MHz. The 
VIR reconstruction using 51 polygons is shown in Figure Hb, where the Voronoi cells have 
also been drawn. Figure |Hi shows the same model but without drawing the Voronoi mesh. 

6.4. Results 

The quality of each reconstruction can be assessed by visual inspection, comparing the 
VIR and MEM model images with the true image. The MEM model looks similar to the 
true image but is noisy. The density of Voronoi generators in the VIR model is greater 
where there is more emission in the true image, approximating the true image with only a 
few polygons. We calculated Xim = J2ii^r°'^ ~ ^i^^'^Yy where 7™°'^ is the intensity at pixel i 
of the model image (MEM or VIR), I^"^^"^ is the intensity at pixel i of the true image, and 
the sum extends over all pixels in the images, xfm gives a measure of how well the model 
fits the true image. It can be seen in Tabled] that the VIR reconstruction has a better Xim 
than MEM, showing that the VIR model is closer to the true image than the MEM model. 



^We choose to display the sky images in a larger field than the domain of free parameters; larger fields 
are required to highlight secondary side-lobes 
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Figures and Hh show the VIR and MEM models residuals. Residual images are the 
dirty map of the residuals of the visibilities, calculated over the optimal model visibilities. It 
can be noted on Figure that the VIR residuals are very good, showing only noise. On the 
other hand, in the MEM residuals (Figure Hh) the object shape can clearly be distinguished 
as well as the CBI's side-lobes. The object seems to be more compact in the model than in 
its MEM residuals; as expected these residuals are convolved with the synthetic beam. 

Restored images are shown in FiguresHf andHi. These images are obtained by convolving 
the models with a Gaussian point spread function (PSF) given by DIFMAP and adding the 
dirty map of the residuals visibilities. On Figures Hf and Hi it can be assessed that VIR 
produces improved restored images relative to MEM. The VIR restored image is similar to 
that expected given the instrumental noise: it approximates the true image convolved with 
a Gaussian PSF plus a uniform noise level. In the MEM restored image, on the other hand, 
the CBI side-lobes can still be distinguished. 

The number of optimization parameters in MEM are 64 x 64 = 4096, while the VIR 
method has only 51 triplets (cell's (x, y) position and intensity) i.e. 153 free parameters. This 
smaller number of parameters causes the Bayesian entropy to be greater than the pixelated 
version, obtaining a smaller value for our merit function L to be minimized. 

2 

Table [1] also shows — — values, where ndata is the number of data points (3620 x 2 in 

2 

our case). A good reconstruction should have a — — value close to 1. It can be seen that 

' ^ "data 

2 

the VIR model gives a value of — — closer to 1 than the MEM reconstruction. 

^ "data 

7. Conclusions 

We have introduced a Bayesian Voronoi image reconstruction (VIR) technique for inter- 
ferometric data where the image is represented by a Voronoi tessellation in place of the usual 
pixelated image. The advantage of Voronoi models is that we can use a smaller number of 
free parameters, as required by the Bayesian analysis of a discretized intensity field. Our 

Table 1. Comparison between MEM and VIR reconstructions. 







"data 


L 


X- 


MEM 


7354.85 


1.016 


12192.6 


0.001608 


VIR 


7221.04 


0.997 


3753.28 


0.001396 
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purpose is not optimal CPU efficiency; we search for the optimal image and model from 
a Bayesian point of view. The free parameters of our model are the Voronoi generators 
positions (xj, yi) and intensities /j. The following points summarize our work: 

• We discretized the intensity field in order to calculate a priori probabilities. We defined 
a quantum intensity value (Xq such that /j = a^N^, where /j is the intensity at cell i 
and Ni the number of quanta in cell i. 

• We calculated the analytical derivatives required by the conjugate gradient and cross 
checked them by finite differences. Because the parameter space in cell generators 
positions is very structured, the positions of the Voronoi generators are difficult to 
change. As initial condition we took a Voronoi tessellation adjusted to an interrupted 
maximum likelihood reconstruction. 

• We simulated a CBl observation over a true image and reconstructed sky images from 
this mock visibility dataset using MEM and VIR. 

• We defined the value of as the estimated noise of the dirty map and searched for 
the optimal number of Voronoi polygons for our example dataset. 

• We finally compared the MEM and VIR models, residuals and restored images. The 
VIR model is closer to our true image than the MEM model. Residuals and restored 
images are also better in VIR than in MEM. We found that VIR model visibihties give 
a better fit to the data than MEM, in the sense that is closer to its expected value. 

We are grateful to Tim Pearson for advice on FFTs and the use of MOCKCBI. G.F.C. 
and S.C. acknowledge support from FONDECYT grant 1060827, and from the Chilean 
Center for Astrophysics FONDAP 15010003. 



A. 



Derivatives 



Our merit function for minimization is 



L 




(Al) 



(A2) 
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So, the derivative of L with respect to any variable x is 

dL _ldx^ _ dS 
dx 2 dx dx 

A.l. Calculation of the Derivatives of 

derivatives with respect to any variable x can be obtain as follows 

9 2^ 



(A3) 



dx 2^ dx \ 2 



fc=i 

k=l ^ ^ ^ ^ 

where its necessary to calculate the model visibilities derivatives with respect to x. 

A. 1.1. Calculation of —^Y~ 
In our Voronoi tessellation representation of the sky image 

Nv 

h 



(A4) 



V{k)^S2lJ A{x)e^'''^^dx, (A5) 



where Ny is the number of polygons, Vj is polygon i and li its intensity. We neglected the 
— — term which is close to 1, but it can easily be included in A(x). After derivation 
and defining fk{x) = A{x)e^'^'^^'^^ we obtain 



dh 



fk{x)d\ (A6) 



Vi 

sin(7r^fcAa;) smjnv.Ay) y ^^^2.^iu,x,+v,y0 (A7) 
AxAi/ J2 Aie2^^("'=^'+'''=^') (A8) 

pixels /eVi 



for small Ax and Ay. 
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A. 1.2. Calculation of ^ — and — 

To evaluate ^ we move the generator Xi an infinitesimal quantity 5x parallel to the x 
axis as in Figure [6l We will calculate 



where AV^ = Vk{xi,- ■ 




Fig. 6. — Voronoi tessellation before and after translating the site Xi by d^- Voronoi gene- 
rators are represented with dots. The solid lines are the polygons before moving afj. The 
dotted lines represent the new polygons after varying Xi. 

It can be seen in Figure [6] that when moving the generator Xj, the only polygons modified 
are Vj and its neighbors. Using this, Eq. IA5I leads to 

AVfc = I'i I fk{x)dx-Ii / fk{x)dx 




(AlO) 



where Vj is the polygon generated by Xi before moving, V- is the same polygon after moving 
Xi, Ji is the set of indices of the polygons that are neighbors to Vj and J'^ is the set of indices 
of the polygons that are neighbors to V^'. 

It can be seen in Figure [6] that 

v. = (v.nv:)u(VAv,nv;), v; = (v.nv;)u(v;\v,nv;), (aii) 
V, = (V, n vj) u (v; n v,), vj = (v, n vj) u (v. n vj), (A12) 
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so, Eq. ETO] is 



A14 



fk{x)dx 



+ V (/; - /,) / fk{S)dx + il'^ - U) / h{x)dx 



(A13) 



In our case the cells' intensities don't depend of the position of the generators, so we obtain 



AVk = V - Ij) ( [ fk{x)dx - [ fk 



[x)dx 



(AM) 



It can be seen in Figure E] that to obtain AVk we must integrate only over the shaded 
regions. For this purpose, for each region between Xi and Xj we will define a coordinate 
system 



s = — cos ajx + sin ajy. t = sin ajx + cos ajy, 



(A15) 



where aj is the angle formed by the —x axis and the edge aij between Xi and xj (see Figure 
[7]). Using this change of coordinates, the integral over the region of interest is 



fk{x,y)dxdy 



fk{s,t)dsdt. 



(A16) 




Xj Sx Xj 



Fig. 7. — Change of coordinates from {x,y) to (s, t). 
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Let Xi = {xi, Hi) be the position of the i cell's generator, Xj = {xj, yj) one of its neighbor, 
and Xi = {xi + 6x, Vi) the site's position after moving it a quantity d^- We define Xq = {xq, yo) 
as the point in the intersection of the segment formed by Xi and x'j and its respective edge 
ttij. The same way, we define Xq = {xQ,yQ) as the point in the intersection of the segment 
formed by x/ and xj and its respective edge a-^. It can be seen on Figure [7] that xq = ^'"^^-^ 
, x[) = Xo + f and y'^^ = yo = 

The edge defined in the new coordinate system by 

s = Sq = —Xq cos aj + ?/o si^ (^j- (Al''') 

In the same way, the edge a'^^ is defined in the original coordinate system by 

y = m{x - x'q) + yo, (A18) 

where 

Xi -\- 5x Xj . , 

m = -. (A19) 

Vj - Vi 

We can define the same line in our new coordinate system as 



s = m't + b', (A20) 



where 



, cos a,- + m sm a, . . . s 

m = ^ ^, (A21) 

sin aj — m cos aj 

b' ^ -^^'o + y^ . (A22) 
sin aj — m cos aj 

This can be approximated to first order in 6^ as 

m ~ 6xMx, (A23) 
b' ^ So + 45,, (A24) 



where 



sin^ ai sin a,- cos a,- / . 

= ^-^ = (A25) 

yj yi Xi Xj 

sin a,- cosoj . . / a ^^^n 

-Da;. = —{sq COS aj + Xi) = —[Socosaj + Xij. (A2o) 



Vj Vi Xi Xj 

The integral in Eq. IA14I using our new coordinate system will be 



J = / fk{x)dx- / Mx)dx (A27) 
' ^(f)e2-^(«^+^s/)(ia;rf?/. (A28) 
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If we use A{x) in the {s,t) coordinate system as a pixelated image, Eq. IA28I will be 

"^fjl f"i't+b' 
1 

I e pixeles de ijl 



pt^l prn't+b' 

X= "^^11 e2'^*("^(^'*)+"f(^'*))M, (A29) 



where tl^^ and t'^j^ are the t coordinate of the beginning and end of the portion of the edge 
aij that intersects pixel /. Developing the previous expression, 



J _ A; / / ^2m{u{-scosaj+tsmaj)+v{ssmaj+tcosaj))^^^^^ (A30) 

Jt^., Jsn 



where we defined 



Yl ^e2"^('''"i+*»^"^2)Ky7^x-, (A31) 



ci = —ucosaj + vsmaj, (A32) 

C2 = usinaj + f cosaj, (A33) 

Kiji = {MJiji + B^) sin (7rc2Atij7) (A34) 

vrc2 



.M^ /sin(7rc2Atij7) \ 
+z— Miji cos (7rc2Ati 7) , 

/ V 7rC2 / 



U^i - (A35) 
^hi = fcl^. (A36) 

In the calculation above we integrated over the fraction of the edge that falls inside pixel 
/ and then summed these integrals over the whole edge Oj. It is also possible to approximate 
the integral of Eq. IA30I as /^i^' g{t)dt = g{tiji)Atiji, which is equivalent to taking the limit 

ijl 

over the integral I of Eq. IA311 lim^ti^^^ol, obtaining 

I = J2 Ai^tviW,,i + 5,)e2-^(*-'^^+^°^i)(5,. (A37) 



I 



We found by direct evaluation that the difference between Eq. IA37I and Eq. IA31I is 
negligible, so, for simplicity, we will use Eq. IA37I Introducing Eq. I A3 71 in Eq. IA14t we 
obtain 

AVk = 5^Y. [(^^ - ^^■) AiAU,i{MAji + S,)e2-(*-'^^+^o^^)] , (A38) 
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so, according to Eq. IA91 the derivative of the k visibihty with respect to the position x of 
polygon i is 



hm 



dxi &x^o 5x 



(A39) 
(A40) 



Similarly, for the derivative with respect to the position y of the i polygon we obtain 



(A41) 



where 



COS a,- sm a, cos a, 
My = ^ = ^ ^ 



3/ 2 



Vj - Vi 



sm a 



3 



COS a 



Vj - Vi 



,Sosmaj -yi). 



(A42) 
(A43) 



A. 2. Calculation of the Derivatives of S 



We defined our entropy as 



S 



In 



n 

\n{N\) - N \n{n) - ^ ln(Ari!^ 



1=1 



n 

In (^r(Ar + l)) - Arin(n) -^In (^r(iVi + l)), 

1=1 



(A44) 
(A45) 

(A46) 



where A^j = ^ is the number of quanta in cell i, N = Ni and T is the Gamma function. 
It can be seen that this function does not depend on the position of the Voronoi generators, 
so 

(A47) 



as _ as _ ^ 

dxi dyi 



Using Weierstrass' definition of the Gamma function 



n=l 



1 + - 

n 
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where 7 is Euler's constant, we can obtain 



d\n{T{z + l)) - . 

l^ = -7 + E^ (A48) 

n=l 



SO, the derivative of S with respect to Jj is 



k=l k=l 



A. 3. Finite Difference Cross Check on the Derivatives 



Numerical calculation of the derivatives by finite differences is not very accurate, in 
particular for the position of the generators. Finite difference derivatives are calculated as 
_ L{x+s)-L(x) ^ ^Yieie 5 is a small displacement of x. In the case of the positions of the 
generators, if 6 is too small, the pixelization of the Voronoi diagram (needed to obtain the 
model visibilities) will not change after the displacement 6. On the other hand, if 6 is too 
big, the generator displacement may cause the function to change abruptly, as explained 
below. That is why we calculated the analytical expression for the derivatives. 

To verify that our derivatives are correctly calculated and programmed, we compared 
our analytical result with a numerical calculation. We created a Voronoi tessellation of 50 
polygons with random positions and intensities and calculated the analytical and numerical 
derivatives using these parameters {xj, yi, h}. For the case of ^ and ^ this numerical cross 
check consists of moving each Voronoi generator a quantity S from -0.1 to 0.1 with an interval 
of 10~^ in units of the total size of the square image. We evaluate the merit function L at 
each position intervals, thus obtaining two sequences {Li}'^^^^ . We then fitted a polynomial 
of order four to the curve defined by each sequence {L^} and calculated the derivative of the 
polynomial at 5 = 0. For the case of ^ we varied the intensity of cell i from —a^ to Cq 
and did the same approximation to a polynomial of order four and calculated its derivative. 
Figure [H] shows this cross check for ^ and Although the derivatives are similar, they are 
not exactly the same for This is caused by the polynomial coarseness fit, as explained 
below. 

Figure [3 shows the curve fit for ^ for three different generators (generator number 37, 
36 and 18 respectively). It can be seen in Figure [9] that the polynomial fit adjusts quite well 
to the function values for polygon number 37, so on Figure [8] both derivatives are the same. 
On the contrary, for polygons number 36 and 18, the fitted polynomial does not resemble the 
function L at 5 = 0, causing a slight difference in their derivatives on Figure El For polygon 
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number 18 the polynomial does not fit the curve at all. This is the main problem of using a 
numerical approximation for the derivatives of {xi}: when two polygons are closer than 6, 
the generator displacement causes the function L to change abruptly (see Figure [TUl) . 

It can be seen that the analytical and numerical derivatives on Figure [8] are almost 
the same. As explained above, differences are produced because there are cases where the 
polynomials do not fit well to the variations of the merit function L (for example, when two 
generators are too close). In an accurate calculation it is necessary to use the analytical 
derivatives. 



B. Fitting a Voronoi Tessellation to an Image 

Once we have a reasonable reconstruction for a pixelated image, we would like to fit a 
Voronoi tessellation to it in order to have a good initial starting point for the CG. This is 
done in an incremental way. 

We start with a mesh consisting in only one polygon. We calculate the error per polygon 

as 

where the sum runs over all the pixels that fall inside polygon i, li is the intensity of that 
polygon and J™ is the intensity of pixel I in the image to be fitted. In each iteration we 
add a new polygon inside the one with the greatest error. The new generator is inserted in 
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5 

Fig. 9. — Examples of polynomial fits, used to determine numerical derivatives of the opti- 
mization function L for a particular generator. Dots represent L vs the polygon displacement 
5va.x and the solid line shows the fourth order polynomial fit to L. A vertical line is drawn at 
5 = 0, where the derivatives were calculated. Top: Generator number 37, with a satisfactory 
polynomial fit. Middle: Generator number 36, the curve is not a good fit at 5 = 0. Bottom: 
2gth generator, the curve is not a good fit because L shows an abrupt variation near 5 = 0, 
which is due to the proximity of another generator. 

the position of the pixel that has the most different intensity value with respect to the mesh 
intensity. 
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Fig. 10. — Translation of a generator close to another. Left: Before moving generator afj, 
polygon i, the darker polygon in the image, is on the left. Right: After moving generator Xi, 
by a displacement of 6, polygon i is on the right of polygon j. When displacing generator xi 
the diagram changes considerably, with a concomitant abrupt variation in L. 
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